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Abstract 

We propose two supersymmetric Standard Models (SMs) with decaying and stable dark matter 
(DM) particles. To explain the SM fermion masses and mixings and have a heavy decay DM 
particle S, we consider the Froggatt-Nielsen mechanism by introducing an anomalous U(l)x gauge 
symmetry. Around the string scale, the U(l)x gauge symmetry is broken down to a Z2 symmetry 
under which S is odd while all the SM particles are even. S obtains a vacuum expectation value 
around the TeV scale, and then it can three-body decay dominantly to the second/third family 
of the SM leptons in Model I and to the first family of the SM leptons in Model II. Choosing a 
benchmark point in the constrained minimal supersymmetric SM with exact R parity, we show that 
the lightest neutralino DM is consistent with the CDMS II experiment. Considering S three-body 
decay and choosing suitable parameters, we show that the PAMELA and Fermi-LAT experiments 
and the PAMELA and ATIC experiments can be explained in Model I and Model II, respectively. 

PACS numbers: 12.60.Jv, 14.70.Pw, 95.35.+d 
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I. INTRODUCTION 



It is well known that supersymmetry privides an elegant solution to gauge hierarchy 
problem in the Standard Model (SM). In the Minimal Supersymmetric Standard Model 
(MSSM), gauge coupling unification can be realized, which give us the important hint of 
Grand Unified Theory (GUT). In addition, in the supersymmetric SMs, we can define a Z 2 
symmetry called R parity under which the SM particles are even while their supersymmetric 
partners are odd. With R parity, we can avoid the dimension-four proton decay problem and 
evade the stringent constraints from the electroweak precision data naturally. Interestingly, 
the lightest supersymmetric particle (LSP) is stable due to the R parity, and then can be 
the dark matter (DM) candidate. For example, the lightest neutralino can be a viable cold 
DM candidate, which can give us the correct relic density as well. 

During the last two years, there were quite a few very interesting DM experiments from 
indirect and direct dections. The ATIC [l] and PPB-BETS 2| collaborations have reported 
the measurements of cosmic ray (CR) electron/positron spectra at energies up to ~ 1 TeV. 
These data show an obvious excess over the expected background in the energy ranges 
~ 300 - 800 GeV and ~ 500 - 800 GeV, respectively. In the mean time, the PAMELA 
collaboration also released their first CR measurements of the positron fraction jsj and the 
p/p ratio Q. Although the p/p ratio is consistent with the astrophysical expectation from the 
interactions between the CR nuclei and interstellar medium, the positron fraction indeed 
shows a significant excess for energies above 10 GeV up to ~ 100 GeV, compared to the 
background predicted by the conventional CR propagation models. Later, the Fermi-LAT 
collaboration has released data on the measurement of the electron/positron spectrum from 
20 GeV to 1 TeV with unprecedented precision and the HESS collaboration has released 
the data on the measurements of electron/positron spectrum from 340 GeV to 700 GeV jsj, 
complementing their earlier measurements from 700 GeV to 5 TeV 0j. For simplicity, we 
will denote the Fermi-LAT collaboration as FERMI collaboration in the following. Although 
the corresponding data from the ATIC and the FERMI/HESS experiments are not fully 
consistent, it was shown that the DM models, where the DM particles annihilate or decay 
dominantly to the SM leptons, can explain these experiments by choosing the suitable DM 
particle mass and the proper final state particles. 

To explain the PAMELA/ ATIC experiments or the PAMELA/FERMI/HESS experi- 
ments from DM annihilations, we know that a large boost factor about 100-1000 is needed. 
However, from astrophysics, the N-body simulation shows that the boost factor from DM 
substructure can never be larger than 10 jqj. To solve this problem, one can consider the 
Sommerfield enhancement js-ll] or Breit-Wigner resonant enhancement 12 1. Alternatively, 
we can also consider the non-thermal dark matter production so that the DM annihilation 
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cross section can be large 13j. In addition, if the DM particle is not absolutely stable and can 
decay dominantly to leptons, we can explain these experiments for the DM lifetime at the or- 
der r ~ 10 25 — 10 27 s 14H18j| . In particular, in the supersymmetric Standard Models, the LSP 



neutralino cannot explain the PAMELA/ATIC experiments or the PAMELA/FERMI/HESS 
experiments unless it can decay due to the suitable i?-parity violation dimension-four op- 
erators. Furthermore, to fit the PAMELA and ATIC data via the Markov Chain Monte 
Carlo (MCMC) technique, one found that the DM mass is about 700 GeV for annihilation 
and 1.4 TeV for decay, and the favored final state is e + e~ 19|. And to fit the PAMELA, 
FERMI, and HESS data, one found that the DM mass is about 2 TeV for annihilation and 
4 TeV for decay, and the favored final states are the combination of fi + fi~ and t + t~ since 
the electron/positron spectra in the FERMI and HESS experiments are softer than these 
in the ATIC and PPB-BETS experiments 19|. Also, the HESS obervation of the Galactic 
center gamma rays gives strong constraint on the annihilation DM scenario while gives much 
weaker constraint on the decay DM scenario. Thus, it favors the decay DM 

Recently, the Cryogenic Dark Matter Search (CDMS) collaboration has observed two 
candidate DM events in the CDMS II experiment 20j. The recoil energies for these two 
events are 12.3 keV and 15.5 keV, respectively, and the data set an upper limit on the DM- 
nucleon elastic-scattering spin independent cross section around 10~ 8 — 10~ 7 pb. Because 
the probability of observing two or more background events is 23%, the CDMS II results 
cannot be a statistically significant evidence for DM interactions, but these two events can 
not be rejected as signal. In particular, the favored DM mass from the CDMS II data is 
about 100 GeV. Later, the CDMS II results have been studied extensively in various DM 
models 21|. Interestingly, the CDMS II experiment can be explained in the supersymmetric 



Standard Models where the LSP neutralino is DM. 

In short, if the PAMELA, ATIC, FERMI, and HESS experiments indeed observed the DM 
annihilations or decays, the corresponding DM particle is heavy around a few TeVs. And 
if the two events observed by the CDMS II experiment are DM signals, the corresponding 
DM particle is light around 100 GeV. Therefore, there may exist at least two DM particles 
in the Nature. In fact, in almost all the previous DM models, the nearly universal implicit 
assumption is that there is one and only one DM particle. However, we cannot prove this 
implicit assumption, and then we cannot ingore the possibility of multicomponent DM 0, 



23|. 

In this paper, we propose two supersymmetric Standard Models with decaying and stable 
DM particles. To avoid the proton decay problem and evade the stringent constraints from 
the electroweak precision data, we assume that R parity is not violated. In our models, we 
require that the LSP neutralino be the stable DM particle with mass around 100 GeV and 
the LSP neutralino-nucleon scattering cross section be about 10~ 8 pb. Thus, we can explain 
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the CDMS II experiment. We also assume that the supersymmetry breaking scale is still 
below 1 TeV, and then we can solve the gauge hierarchy problem without fine-tuning. To 
explain the PAMELA, ATIC, and FERMI experiments, we introduce a DM particle S with 
mass around a few TeVs. To produce the SM fermion masess and mixings and have the 
heavy decay DM particle S, we consider the Froggatt-Nielsen (FN) mechanism 24J. We 
introduce an anomalous U(l)x gauge symmetry whose anomaly is cancelled by the Green- 
Schwarz mechanism 25]. Especially, the U(l)x gauge symmetry is broken down to the Z 2 
symmetry around the string scale under which S is odd while all the SM particles are even. 
Thus, S can be a DM particle. Similar to the discussions in Refs. 
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2911 , the SM fermion 



masses and mixings can be generated as well. With a pair of heavy vector-like particles 
that are SM singlets and have U(l)y charges ±1 respectively, we obtain that the leading 
Yukawa coupling terms between S and the SM particles are S 2 HdL i E'j where Hd is the Higgs 
field, and Lj and E c - are the 2-th family of the SM lepton doublet and j'-th family of the 
right-handed charged lepton, respectively. In Model I, we have (i, j) = (2, 2) and (3, 2), 
while in Model II, we have (i, j) = (1, 1). After 5* obtains a vacuum expectation value 
(VEV) at the TeV scale, S can decay dominantly to the second/third family of the SM 
leptons via dimension-six operators (three-body decay) in Model I, and to the first family 
of the SM leptons in Model II. To realize our idea, we present a benchmark point from the 
constrained minimal supersymmetric Stanard Model (CMSSM). The lightest neutralino N® 
contributes to part of the whole DM relic density, i.e., VL^oh 2 « 0.08. The LSP neutralino- 
nucleon elastic-scattering spin independent cross section is about 5 x 10~ 9 pb. Thus, this 
benchmark point is consistent with the CDMS II results. In addition, for the S lifetime 
about r ~ 10 25 - 10 27 s, we can explain the PAMELA, FERMI and CDMS II experiments 
in Model I with S mass 3 TeV, and explain the PAMELA, ATIC and CDMS II experiments 
in Model II with S mass 1.8 TeV. 

The paper is organized as follows. In Section II, we explain the SM fermion masses and 
mixings as well as the leading Yukawa terms S 2 LiHdEj via the FN mechanism. In Section 
III, we present a benchmark point in the CMSSM parameter space, and we discuss the 
DM particle S three-body cascade decay. We fit the PAMELA and FERMI data and the 
PAMELA and ATIC data in Section IV. Section V is our discussion and conclusions. We 
present more technical details in Appendices A, B, and C. 



II. TWO SUPERSYMMETRIC STANDARD MODELS 

In this Section, we will present two supersymmetric Standard Models with decaying and 
stable DM particles. To have a decay DM particle, we introduce a SM singlet field S 
and a new Z 2 symmetry. Under this Z 2 symmetry, S is odd while all the SM particles 
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are even. Thus, S can be a DM particle. Because the global discrete symmetry can be 
broken by the quantum gravity effects, it is natural to have the remnant discrete symmetry 
from an extra U(l)' gauge symmetry |30|, |3l|. In particular, in the U(l) '-extended Minimal 
Supersymmetric Standard Model (UMSSM) 32] where we can solve the /i problem and 
avoid proton decay, we indeed have a new DM candidate if the U{1)' is broken down to a 
Z 2 symmetry j^, 33 1 . 

In our models, this Z 2 symmetry is the residual symmetry after the anomalous U(l)x 
gauge symmetry breaking, i.e., Z2 is a subgroup of U(l)x- After supersymmetry is broken 
at the TeV scale, the supersymmetry breaking soft masses may make m| negative and then 
the scalar component of S obtains a VEV around the TeV scale. We assume that the 
lightest state in the supermultiplet S is its scalar component, which is still denoted by S 
in this paper. In short, S is a DM particle which can decay after Z2 symmetry breaking. 
In the concrete model building, we may need a singlet sector with additional particles and 
interactions for two kinds of reasons: (1) We have to stabilize the VEV of S at the TeV 
scale; (2) We need to generate the correct S relic density from thermal productions by 
considering the additional U (1)' gauge symmetry at the TeV scale 23L or from non-thermal 
productions by considering the cosmic string at intermediate scale 13( or the extra vector- 
like particles at the intermediate scale 16j. Similar to the Ref. 16], we concentrate on the 
DM phenomenology in this paper, and then the detailed discussions of the singlet sector and 
S relic density are out of the scope of our current paper. 

To explain the PAMELA, ATIC, FERMI and HESS experiments, S must decay dom- 
inantly to the SM leptons. To select the suitable final state particles, we consider the 
Froggatt-Nielsen mechanism 24| , which is an elegant way to explain the SM fermion masses 
and mixings. We consider an anomalus U(l) y gauge symmetry, whose anomaly can be 
cancelled via the Green-Schwarz mechanism |25|. In the string model building, we indeed 



have at least one anomalous U(l)x gauge symmetry 29]. 

To break the U(l)x gauge symmetry, we introduce a flavon field A with U(l)x charge 
— 1. Because supersymmetry must be preserved close to the string scale, A can acquire a 
VEV so that the U(l)x D-flatness can be realized. It was shown 27| that 

(A) 



0.171 < e 



Mi 



< 0.221 



pi 



(1) 



where Mpi is the reduced Planck scale. Interestingly, e is about the size of the Cabibbo angle. 
Let us explain our convention. We denote the SM left-handed quark doublets, right-handed 
up-type quarks, right-handed down-type quarks, left-handed lepton doublets, right-handed 
neutrinos, and right-handed charged leptons as Qi, Uf, Df, Li, N?, and Ef, respectively. 
Also, there is one pair of Higgs doublets H u and Hd in the supersymmetric Standard Models. 
Moreover, we introduce a pair of vector-like particles E' and E' whose quantum numbers 
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under SU(3)c x SU(2)l x U(1)y are (1, 1, —1) and (1, 1, 1), respectively. In addition, the 
U(l)x charges for all the particles in our models are denoted by appropriate subscripts, for 
example, for a generic particle <ft, its U(l)x charge is X^. In our models, we choose the 
U(l)x charges for the SM fermions and Higgs fields as follows 
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The superpotential for the SM fermion Yukawa couplings is 
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(3) 



m p1 ; V m pi, 

where i, j = 1, 2, 3, and all the coefficients are assume to be order one in the above 
superpotential. Thus, we obtain the SM fermion Yukawa coupling matrixes Y u , Yd, Y e and 
Y u respectively for up-type quarks, down-type quarks, charged leptons and active neutrinos 



(4) 



where T is tr ansp ose. We can show that the observed SM fermion masses and mixings can 
be generated (27|, bd], and tan/3 = (H®)/(H d )° is 25. 

In the Model I, we would like to explain the PAMELA and FERMI experiments. We 
choose the following U{l)x charges for S, E' and E' 





/ £ 8 £ 5 £ 3 ^ 




/ e 5 e 4 e 4 ^ 




/e 2 


6 e\ 


Y u ~ 


e 7 e 4 e 2 


, Y d ~ Yj ~ 


e 4 e 3 e 3 




e 


1 1 




^ e 5 e 2 1 J 




e 2 e 1 e 1 ) 




V c 


1 1/ 



A 5 — 3/2, A S / 



A 7 



5. 



(5) 



Note that A has U (l)x charge —1, the U(l)x charges for the SM particles are integers, while 
the U(l)x charge for S is half integer, thus, the U(l)x gauge symmetry is broken down to a 
Z2 symmetry after A obtains a VEV. In particular, under this Z2 symmetry, only S is odd 
while all the other SM particles are even. Then the leading Yukawa coupling terms between 
S and the SM particles in the superpotential are 

5 
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E C E' + M V E'E' 



(6) 



where k = 2, 3, and My is the vector-like particle mass around the 10 13 GeV. In particu- 
lar, the superpotential mass terms EfE' and A n E^E' are forbidden due to the holomorpic 
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property of superpotential. Because E' and E' are heavy, we should integrate them out 
below their mass scale. From the above superpotential, we obtain from equations of motion 
(EOMs) for E' and E 1 

E> = - , E> = - S2E * . (7) 

My ' MyMpi V ' 

Using the above EOMs for E' and E', we obtain that the superpotential in Eq. ([6]) becomes 

W ^ ~ e5 ( 1 H d L k E c 2 . (8) 

\M m M v ) 2 K 1 

To simplify the discussions, we assume that the coefficient of S^H^LiE^ * s about three 
times larger than that of S^H^L^E^- Thus, we will concentrate on the term S 2 H ( iL 2 E2 in 
the following discussions. 

In the Model II, we want to explain the PAMELA and ATIC experiments. We choose 
the U{l)x charges for S, E' and E' as follows 

X 5 = 1/2, X E , = -A, X^=A. (9) 

Similar to the Model I, the U(l)x gauge symmetry is broken down to a Z 2 symmetry after 
A obtains a VEV. Under this Z 2 symmetry, only S is odd while all the other SM particles 
are even. Then the leading Yukawa coupling terms between S and the SM particles in the 
superpotential are 

w d HdLiW + (£) E;e + MvWe ■ (io) 

Similar to the Model I, integrating out the vector-like particles E' and E', we obtain 

wd -« 4 (s!k) w - (u) 

In this paper, we will define M 2 = M^\Myje n where n is equal to 5 in Model I and 4 in 
Model II. With My around 10 13 GeV, we obtain that M* is around 10 17 GeV. To get My 
around 10 13 GeV, we can introduce additional global U(l)" symmetry and a SM singlet field 
S' which breaks the U(l)" symmetry. The U(l)" charges for Li, E![, E', S, and S' are 1, — 1, 
— 1, 1/2 and 1, while the U(l)" charges for all the other particles are zero. Therefore, all 
the previous terms in the superpotential including the SM fermion Yukawa coupling terms 
are invariant except the term MyE'E'. The vector-like mass term for E' and E' can be 
generated by the following superpotential term after the U(l)" symmetry breaking 

W D S'WE' . (12) 
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Assuming that S' acquires a VEV around 10 13 GeV, we obtain that My is around 10 13 GeV. 

By the way, the other Yukawa coupling terms among S, the SM fermions and Higgs fields 
might have the prefactor ( -j^r) ijj^j (•••) (n > 1 or m > 1). Thus, these Yukawa cou- 
pling terms are supressed by at least the Planck scale square and then are negligible since S"s 
VEV is around the TeV scale. In this paper, we only consider the minimal Kahler potential. 
In general, the dimension-six operators of the form SS'cjx/y /M^^ cannot be forbidden by any 



hese dimension-six operators 
1^ | . However, their contribu- 



symmetry where <fi denotes the SM fermions and Higgs fields. ' 
can induce DM two-body decay through derivative couplings 
tions to the CR are also ignorable since these operators are suppressed by the Planck scale 
square as well. 

III. DECAY AND STABLE DARK MATTERS 

A. The CMSSM Benchmark Point for CDMS II Experiment 

The CMSSM with the LSP neutralino as DM has been studied extensively before. Typi- 
cally, the LSP neutralino with mass about tens of GeV is in tension with the lightest CP-even 
Higgs boson mass whose low bound from the LEP is 114 GeV. Especially, we require that 
the LSP neutralino relic density be small about half of the total DM relic density, i.e. 
Qfioh 2 ~ 0.06. In this paper, we do not scan all the viable parameter space. We only con- 
sider a CMSSM benchmark point that satisfies all the constraints. We choose the following 
five free parameters at the GUT scale 

m = 310GeV, m 1/2 = 250GeV, A = -1040, tan/3 = 30, signQi) = +1. (13) 



As expected, this benchmark point is in the coannihilation region 34j with small mass 
difference between the light stau T\ and LSP neutralino, i.e., — m^o < 8 GeV. Such 
benchmark point is not interesting previously since the LSP neutralino relic density is smaller 



than the observed whole DM relic density, 
two DM particles. With MicrOMEGAs 2.0 



dowever, it is fine in our models since we have 



351 ] . we obtain the LSP neutralino relic density 



%o/i 2 ^0.08. (14) 

Moreover, the mass of the LSP neutralino is about 101.6 GeV. And the spin-independent 
cross section between the LSP neutralino and nucleon is about &si ~ 5 x 10~ 9 pb. Although 
this cross section is a little bit small, we can still explain the CDMS II experiment due to 
the uncertainties of the QCD effects in the calculations. 

It is necessary to address the mass spectrum of this benchmark point in details for the 
following calculations. The lightest neutralino N% is bino-like, while the heavy neutralinos 



S 



^3,4 are Higgsino-like. Since tan f3 is large, in the charged Higgs boson system, is the 
major component of the Goldstone boson, with the ratio to (Hj)* to be tan 2 /3 : 1. So we 
can treat Hj as the charged Higgs boson C~ approximately, i.e., H~[ ~ C~ . The CP-odd 
neutral Higgs and charged Higgs mass matries are diagonlized by the same matrix. Thus, the 
imaginary part of H® gives the main component of the Goldstone boson, while ^4° is mainly 
the imaginary part of H®, i.e., A « Im(H^). Same conclusion applies to the heavy CP-even 
Higgs boson H°, i.e., H° « Re(H%) because H° is equal to cosaRe(if^) + sin aRe(H®) with 
a ~ —0.04. As for the chargino system, the large /i term implies that the charged Higgsino 
is the major component of heavy chargino C^- Because we need to use the neutralino, 
chargino and Higgs mixing matrices in the following discussions, we calculate the relevant 
masses and mixings at low energy by SuSpect 36]. For the neutralino system, the four 
particles (—iB,—iW 3 ,H%,H®) are transformed into the mass eigenstates (iV°, iV° , iVg, iVj) 
by the unitary matrix Z 

( -iB \ ( 1.00 -0.02 0.08 -0.02 \ 
0.03 0.99 -0.16 0.06 
0.04 -0.07 -0.70 -0.71 
0.07 -0.15 -0.69 0.70 



-iB \ 
-iW 3 

V H2 J 



\ 



N° 3 



I 



(15) 



And the neutralino mass eigenvalues are m 



1*1,2,3,. 



101.6, 196.2, 571.0, 577.0 GeV, respec- 
tively. The particles {W 3 ~ , H~[) and (Re(H®), Re(H®)) can be written respectively in terms 
of their mass eigenstates as follows 



Wi 

h: 



-0.98 0.22 
0.22 0.98 




( MH° d ) 



1.00 0.04 
-0.04 1.00 




(16) 



196.1, 578.5 GeV, and m h o )H o = 117.90, 510.81 GeV. The lm{H%) and lm{H°) 
are not given explicitly here, and the mass of the CP-odd Higgs field A is 512.49 GeV. 



with Vflr,- 

°1,2 



B. DM S Three-Body Decays 

In this paper, we will concentrate on the calculations in Model I. For simplicity, we will 
not explain the calculations in Model II since the calculations are not only similar but also 
simpler. To analyze the DM S primary decays, we write the operator in the components 
explicitly 

C^SH d LE c +c.c. = C^SH° d EE c - C ll SH^uE c + c.c. 

+ C^H^lJi^ + H^Jl L fi R + H%Jil'Pr + #2#£/4) 

- C^H^v,^ + H^fi* R + H-fiA + H d "A), (17) 
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where C M is equal to (S)/M 2 times the order one coefficient. If these states are transferred 
to the mass eigenstates, further mixing factor must be included. After H® obtains a VEV, 
we shall have two-body and three-body S decays. Comparing to the two-body S decays, 
the three-body S decays have an extra factor (Md/ (H®)) 2 / (96tt 2 ) where mo is the mass of 
the heavy DM particle S. In our models, we choose Md = 3 TeV and tan/3 = 30, and then 
we get (Md / (H®)) 2 / (96tt 2 ) ~ 282, Therefore, we only consider three-body S decays in this 
paper. 

The heavy decay DM particle is the scalar component of S, whose three-body on-shell 
decays are given by the second line of Eq. (fT7|) 

S ->■ N 3A [iJl R , Cx^y^L, C lt2 ^ L , A°(H°)nn, C~v^ R . (18) 

With the Eq. flB9j) in the Appendix B, we obtain that all these nine primary decay channels 
almost have the democratic rates except the mixing factor for each channel. Anyway, the 
total decay rate is given by 



1 (S) 2 m 2 D 
total 128tt 3 Mf 



CjRA xm D , (19) 



where I denotes the 7-th decay channel in Eq. ( [18]) . In our paper, the DM particle S is so 
heavy that the phase space suppressing factor Ri is almost process independent and is a 
constant around 1/6 (see Eq. flB9j) ). Then, the scale of the DM S lifetime is estimated to be 

The random coefficient of the operator has been set to be 1. Notice that except the last 
factor, Eq. ( 120]) is the so-called r e // which is the inverse of F e ff defined in Eq. flB9j) . 



IV. COSMIC RAY ANOMALIES 



Although the very large astrophysical uncertainties do exist in the calculations of cosmic 
ray, we do not want to scan all the viable parameter space. We only want to demonstrate 
that with appropriate parameters, we can explain the PAMELA and FERMI experiments 
in Model I and explain the PAMELA and ATIC experiments in Model II. 



A. Cosmic e Excess 



In the Micky Way dark halo, the DM S decays into the SM particles, which propagate to 
the solar system. The propogation of charged particles is described by diffusion equation. 
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For instance, the diffusion equation for positron is 

^ - V • (K(x, E)V^) - A ( b ( E )i>) = q(x, E), (21) 

where ip(x,E) denotes the positron number density per unit energy. Diffusion coefficient 
K = KqE 5 is space independent. The third term describes the energy loss of positron 
through the synchrotron radiation and inverse Compton scattering, with loss ratio b(E) = 
E 2 /te where te = 10 16 s. q(x, E) is the positron source term, describing the positron number 
density injected at x, per unit time and energy. For simplicity, it can be expressed as the 
product of an astrophysics factor and a particle physics factor 

yr s dE J m D 

where we introduce another index F to distinguish the different states from the 7-th decay 
channel. Moreover, the Milky Way DM density profile p(x) in the spherical coordinates is 
generically written as 

^=P°{t) { 1 + (r/r s r ) ' (23) 

where r Q = 8.5 kpc is the distance from the solar system to the Milky Way center, and the 
DM profile density in the solar system is set to be p Q = 0.3GeVcm -3 . According to the 
N— body simulation, it has three popular choices parameterized by (a,/3,j). We use the 
NFW profile with (1,3, 1) 37|, and choose the DM central core r s = 20 kpc. By the way, 
the injected source term is linearly proportional to the decay DM relic density. Therefore, 
the fraction of the DM S relic density rp>M to the whole DM relic density should be taken 
into account. For simplicity, we absorb it into the redefinition of S lifetime by t$ — > t^mTs 
without disturbing any other terms. 

The normalized energy spectrum distribution function dNf F /dE = (2/mD)dNf F /dx 
describes the positrons from the specific DM decay intermediate state (I, F), with (I, F) = 
JiL, A ..., and so on. And Bj p = Bj is the corresponding branch ratio. In this paper, 
summing over all the channels (see Appendix, with a little bit changed notation), we obtain 
the total fragmentation function 

(CI f R,,f) x — — » -j?, (24) 



dx l total ' Rl,F dXj p dx e 



where dN/dx^p and l/Rj t pdN/dx^ ,F respectively give the fragmentation functions in the 
rest frame of the DM S and particle (I, F). The former is analytically calculable, while the 



later is obtained through PYTHIA 38[ simulation (except for muon decay, we shall use the 
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analytical formula). The total spectrum is their convolution, whose explicit form can be 
found in Eq. (IB 1 2 [) . By the way, when we put Eq. 024)) into Eq. f[2"2"j) . the t$ factor in the 
front of source term in Eq. (122]) is cancelled by factor T to tai in Eq. ( |24"|) . So, r e // or T e ff is 
still the free parameter in the calculations of positron propagation. 

Eq. (1211) can be approximately solved analytically. The final positron flux observed at the 
Earth $ e +(r Q , E) is factorized into the convolution of 'H e (E s ,E) which encodes the whole 
astrophysical information and the inject spectrum dNg/dE s 



e 



$ e+ (r ,£) = '^-Kf 2 j E n e (E s ,E)dE sl ^ , (25) 

where /3 e + is the velocity of the observed positron and k = p & / {r e ff mo)- The concrete 
method to calculate the halo functions for electron l-L e (E s ,E) and for anti-proton 7ip(E) 
can be found in Ref. |g[ ■ 

To compare the data, we should consider the cosmic background. It is believed that 
the astrophysical source is mainly due to supernova explosions (primary e~) and the 
interactions between the CR nuclei and light atoms in interstellar medium (secondary e^). 



They are usually parameterized in the following form 40] 



^bk,pri 0.16E 11 ^bk,sec _ 0.7E - 7 



1 + HE - 9 + 3.2E 2 - 15 ' 1 + 110E 1 - 5 + 600£ 2 - 9 + 580£ 4 - 2 ' 

^bk,sec t.a±j ('26") 

e+ ~l + 650£ 2 - 3 + 1500£ 4 - 2 ' 1 ' 

where the unit is GeV~ 1 cm~ 2 s~ 1 sr~ 1 . (j> ^bk,sec arg ^ e background fluxes 

for the primary electron, secondary electron, and secondary positron, respectively. Then 
the positron fraction observed by the PAMELA experiment and the normalized total elec- 
tron/positron fluxes observed by the FERMI and ATIC experiments are respectively given 
by 

<f> , Qbk,sec qDM 

$ + + $ K $ bk 'P ri + fobk,sec ^bk,sec ^ DM $DA/ ' 

e e+ e~ e+ e~ 

E 3 e ($ e+ + $ e -) = E 3 e (K^ pri + $^ sec + $^' sec + $f + A/ + $f A/ ) , (27) 

where k < 1 is a parameter which includes the uncertainties of the primary background 
electron production. 

The PAMELA positron fraction data and FERMI electron/positron data can be fitted 
very well by choosing proper T e ff and k. The primary DM S three-body decays into muons, 
which subsequently decay to electron/positron, give the hard positron energy spectrum 
that is necessary to explain the steep rise in positron fraction observed by the PAMELA 



experiment 17J. As for the FERMI experiment, its data show a rather flat spectrum up to 
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10 20 50 100 200 500 1000 1° 20 50 100 200 500 1° 00 



FIG. 1: The comparison between the contributions to total e fluxes from \i in Model I (or e in 
Model II) and all the other particle cascade decays. The red line denotes the contribution from \i 
(or e), while the dashed line denotes the contributions from the other particles. The left figure is 
for Model II with = 1.8 TeV, and the right figure for Model I with mu = 3.0 TeV. 

about 1 TeV. Moreover, there are a faint minimum and a faint peak at about 100 GeV and 
400 GeV, respectively, which can be considered as the subtle fine structure. Interestingly, 
we can fit such structure as well. The possible reason is that: at low energy (< 50 GeV), 
the DM S decay contributions to the observed electron/positron spectrum are dominated 
by these from A°,H°... channels; while at high energy, they are dominated by these from 
cascade muon decays. The muon decays produce a relative harder spectrum, which is peaked 
around 400 GeV. Because the peak varies with the DM particle mass, this explains why we 
choose mo — 3 TeV as a typical fitting value in Model I. In Fig. (TJ we compare these two 
contributions in Model I with mo = 3.0 TeV and in Model II with mo = 1-8 TeV. In Fig. |2j 
we present the fitting for the FERMI data in Model I with r e// = 0.72 x 10 26 s, m D = 3.0 
TeV, and k = 0.65, and the fitting for the ATIC data in Model II with r eff = 0.52 x 10 26 s, 
van = 1-8 TeV and k = 0.67. In short, these fittings are pretty good especially at relative 
low energy. 

In addition, from Fig. [3J we obtain that our fitting for the PAMELA positron fraction 
data in Model II looks fine since the positron spectrum from the DM direct three-body 
decay is relatively harder. However, we can not fit the PAMELA positron fraction data in 
Model I very well in the low energy region (< 20 GeV). The reason is that the low energy 
positron data are dominated by the sparticle and Higgs boson cascade decays, which produce 
large number of soft electron/positron. Therefore, the positron fraction can exceed the 
background value. Interestingly, the main character at the high energy region is reproduced 
pretty well. Especially, the steep rising is obvious due to the harder positrons from muon 
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FIG. 2: The FERMI data fitting in Model I and the ATIC data fitting in Model II. 

decay. Moreover, it predicts a continuously going up trend until arriving at the peak around 
0.5 TeV, which can be tested by the upcoming PAMELA experiment. 

B. Anti-Proton Fluxes 

Bcause Higgs bosons and charginos can couple to the SM quarks, their decays may 
produce anti-proton excess as well. However, the anti-proton excess is not observed by the 
PAMELA experiment. The propagation of anti-proton is similar to that of positron, but it 
is simpler by the virtue of no energy loss process from the inverse Compton scattering 

^ - V ■ (K(x, E)V4j) + V ■ (V c (x)4j) = q(x, E) - 2h5(z)T ann ^ . (28) 

However, the Milky Way's galactic wind V c (x), which is usually assumed along the z axial 
direction and then is reduced to V c sign(z)k, induces the drift of anti-protons during prop- 
agation, as reflected in the third term. The last term represents the annihilations between 
the anti-protons and interstellar protons in the galactic plane, whose thickness is about 
h ~ 0.1 kpc <C L, the half-thickness of cylinder diffusion region of charged particles. The 
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FIG. 3: PAMELA data fitting for positron fraction in Model I and Model II. 



annihilation rate is given by T ann = (nn + 4?' 3 nHe)&pp n Vp with 41] 

661(1 + 0.0115T- - 774 - 0.984T 00151 )mbarn, for T < 15.5 GeV 



a 



pp 



{ 



36T mbarn, 



for T > 15.5GeV 



(29) 



where T = E — m p is the proton kinetic energy. We can use Kq, 5, V c and L to classify the 
astrophysical models that are compatible with the B/C ratio into three types where K 
and S appear in the diffusion coefficient K. And each type is characterized by the flux of 
anti-protons. To suppress the anti-proton flux, we choose the parameters which can produce 
the minimal anti-protons: 5 = 0.85, Kq = 0.0016 kpc 2 Myr _1 , L=l kpc, and V c = 13.5 km 



Similarly, the anti-proton flux at the heliosphere boundary is formally solved to be 



*!(£) 



1 ^ x d S,H„(E) 



total 



4tt dE 



(30) 



where 



dN 



dN P _r eff y , C 2 Rir) x 1 dN g 
dx V total VTf F Rlf dxr F dxl' F 



(31) 



For details, please see Appendix B. 
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It is much simpler to obtain the flux for anti-proton than that of positron. The point is 
that the astrophysical information is totally factorized out in 7ip(E), and no further convolu- 
tion operation is needed. The predicted anti-proton fluxes comparing to the background flux 
in Model I and Model II are plotted in Fig. HJ where the suitable astrophysical parameters 
have been chosen. The background anti-protons mainly arise from the collisions between 
the primary CR protons (produced also by supernova) and the interstellar hydrogen gas (for 
details, please see Ref. [411]) . Taking proper astrophysics parameters in the practical calcu- 
lations, we find that the anti-proton excess is definitely excluded below 50 GeV (Sixteen of 
the PAMELA seventeen points on p/p are in this energy region j^].). In Model I for the 
PAMELA and FERMI data fitting, the DM particle S is very heavy about 3 TeV, and then 
the anti-proton flux is suppressed in all energy region below about 130 GeV. However, in 
Model II for the PAMELA and ATIC data fitting, the DM particle S is about 1.8 TeV, and 
then the anti-proton flux begins to exceed background above about 60 GeV. Especially, the 
excess is more significant as the energy increasing, which is interesting since it can be tested 
in the future. 



0.0005 



5e-005 



5e-006 



FERMI+PAMELA fit 
ATIC+PAMELA fit 
Background 



10 



50 100 



FIG. 4: PAMELA data fitting for p/p ratio in Model I and Model II. The anti-proton fluxes $%(E) 
is in the unit GeV _1 m _2 s _1 sr _1 . 
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V. DISCUSSION AND CONCLUSIONS 



We proposed two supersymmetric Standard Models with decaying and stable DM parti- 
cles. To explain the SM fermion masses and mixings, we considered an anomalous U(l)x 
gauge symmetry whose anomaly is cancelled by the Green-Schwarz mechanism. Around the 
string scale, the U(l)x gauge symmetry is broken down to the Z 2 symmetry under which 
only S is odd. Thus, S is stable and can be a DM candidate. After S obtains a VEV around 
the TeV scale, the Z 2 symmetry is broken and then S can decay. In our models, on the one 
hand, the LSP neutralino has mass 101.6 GeV, and its relic density is about Qj^oh 2 ~ 0.08. 
Because the lightest neutralino-nucleon cross section is 5 x 10~ 9 pb, the CDMS II results 
can be explained. On the other hand, S is a decay DM particle that can three-body cascade 
decay into the MSSM particles. With suitable U(l)x charges, S decays dominated into the 
second family of the SM leptons in Model I and into the first family of the SM leptons in 
Model II. In Model I, if the mass of the DM particle S is 3 TeV and S has effective lifetime 
about r e ff ~ 0.72 x 10 26 s (the real lifetime is about 2.2 x 10 26 s), we are able to explain the 
PAMELA/FERMI experiments simultaneously. In Model II, if the mass of the DM particle 
S is 1.8 TeV and S has effective lifetime about r e ff m 0.52 x 10 26 s, we can explain the 
PAMELA/ATIC experiments as well. In addition, taking the proper astrophysics parame- 
ters that produce the minimal anti-protons, we have shown that our models are consistent 
with p/p measurement in the PAMELA experiment. 

Finally, although the accompanied gamma-ray and neutrino fluxes are not discussed here, 
they deserve further study. The primary hard neutrinos, as well as the soft gammas and 
neutrinos during the sparticle and Higgs fragmentations are produced. They may provide 
the spectrum properties through the inverse Compton scattering 42] that can be tested by 
the ongoing FERMI experiment. In addition, the complete and systematic studies of the 
multi-body decays and multi-contributions are interesting since we do need them in the DM 
model building. Furthermore, it is interesting to study the CMSSM parameter space where 
the LSP neutralino relic density is smaller than the observed total DM relic desity since we 
may have multicomponent DM. 
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Appendix A: The DM S Three-body Decay 

In this Appendix, we develop some useful approximation formulae for the three-body DM 
S decay. First, let us explain the convention. The decay width is calculated in the rest frame 
of the heavy parent particle S (Fermionic or scalar component depends on the choice.). In 
a decay channel, we always have the order of the final state masses as m 3 > m 2 > vn\. 

In this case, in general m% is the mass of a slepton or Higgs bonson, with m^/mp < 0.1. 
And m 2 denotes the mass of Higgsino, m 2 /mD < 0.05. Finally, the lightest particle is 
either muon or neutrino, whose mass mi = m^ v < 10~ 4 m£> is ignorable. To simplify the 
calculations we set mi = 0. 

First, we consider a general three-body decay process described by a heavy complex scalar 
S with mass M that decays into two fermions ipi^ plus a scalar 3 . The Feynman diagram 
is given in Fig. [51 It is not difficult to get the corresponding amplitude square with the final 



s 




FIG. 5: The DM S three-body decays into the Higgs bosons and leptons, or the sleptons, Higgsinos 
as well as leptons. 

state spin summed over 

\iM\ 2 =AC 2 Cj{p 1 -p 2 + m 1 m 2 ) , (Al) 

where C is the coupling constant with mass dimension —1, and Cj is the possible mixing 
factor arising from converting the interaction state into the mass eigenstate. The mass mi 
is still kept in the momentum, and then will be set to zero finally. 

Appendix B: Fragmentation Function 

Fragmentation function dN/dxi is analytically calculable for three-body decay Y — > 
a + b + c. For the isotropic decay it is very convenient {^J to write the phase space integral 
according to the scaled variables xj = 2Pi ■ Y j VY^, where the invariant mass is defined as 
Y 2 = (a + b + c) 2 . Here, the particles and their four-momenta are labeled with the same 
symbols. In the rest frame of Y, we have Y = my and xi = lEijmy. The three-body phase 
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space is reduced to the two-body phase space by pairing a with b and further integrating 
over c 

Y 2 

d 3 (PS, Y -> abc) = —-^ dx a dx b . (Bl) 

Then we find that the spectrum of the observed final state, for instance a, is obtained by 
integrating over x^. Kinematics constrains x a to be in the region 

2e a ^ x a < 1 + (el - e\ -e\- 2e b e c ) = R™, (B2) 

where e 2 = m 2 a /Y 2 . If all the final state masses are ignorable, the region is between and 
1. For fixed x a , x b is constrained within the region [R^, R^] 

Rt =\{l ~ *a + el)- 1 ^ - x a )(l + e 2 a + e 2 -e 2 c - x a ) ± ^/xJ^Ae 2 ^ (l + d ~ 4, 4)} 



±((2-x a )±x a ) 



(B3) 



where the approximation is valid for small mass limit e a ^, c *C 1, and the triangle function is 
defined as X(x, y, z) = x 2 + y 2 + z 2 — 2xy — 2xz — 2yz. 

The general process is S — > 01^203 with amplitude given in Eq. (lAlj) . In terms of x and 
e, we can write the amplitude as follows 

\M\ 2 =2{CC X ) 2 Y 2 {1 + 4- (e 2 + ei ) 2 - x 3 ) , (B4) 

which is symmetric for the subscripts 1 and 2. The distribution function of the fermion ipi 
is obtained by paring ipi and 03, i.e., setting a = ipi, b = 03, c = -02, 

dN 2(CC X ) 2 Y 4 h . f R t^ J , 2 . N2 . 

x AM ob 3 (1 + e 2 . - (e 2 + Cl ) 2 - x 3 ) 



dx\ 256ir 3 mD 

-J\f 



2(CC X ) 2 Y 4 



1 + el - (e 2 + e x ) 2 - J (i? 3 + (^i) + #3~M 



2567r 3 m/) 

x (tf+OrO-i^i)), (B5) 

where A/" = 1/r is the normalization factor. In the massless limit, it is approximated to 
be a simple function x\. So the fermion spectrum is very hard. The distribution of ip2 is 
got simply by replacing the subscript 1 by 2 in the above formula. Similarly, we get the 
distribution function of scalar 03 by choosing a = 03, b = ipi, and c = ip2 



dN_ _ 2(CCx) 2 Y 
dx 3 256n 3 mD 



/ dx x (1 + el - (e 2 + ei) 2 - x 3 ) 

JR7(xi) 



2 f^Z D M (1 + 4 - {e2 + ei)2 - xs) ( ^ (xs) - • (B6) 
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Its massless limit is simplified to be dN/dx% — > (1 — £3)0:3, so the scalar spectrum is softer 
than the fermion spectrum. Integrating over z 3 , we get the decay rate 

r ' "T^i— x ( c ^) 2y4 x Ris > ( B7 ) 

where 

7^ =1 (1 + 9e 2 3 + 12e* log e 3 + 12e 2 e 3 - 3e*) . (B8) 
o 

In this approximation, the terms involving €\ is negligible here. Furthermore, if e 2 is also 
small enough, the last two terms can be discarded. 

For convenience, we define the DM effective decay rate for all channels as the typical 
scale of its decay 

Then we have Tj = T e ffCjRi f- And the normalized distribution function is rewritten as a 
clear form dN/dx = 1/Rj t p x (dN/dx). 

1. General Cascade Decay Spectra 

In our models, e ± signals come from both the DM cascade decay and direct decay, for 
example, 

S — > Ni^fjfi, /i — > evv, Jl — > N\e. (BIO) 

The positron energy spectrum in rest frame of S is determined to be the convolution of two 
distribution functions. Physically speaking, the fra gme ntation functions of the (N)MSSM 



particles (7, F) can be extracted out from PYTHIA 38 



dN IF 2El> F 

x 1 / = (Bll) 



dX e ™>I,F 

which is calculated in the rest frame of (7, F) particle. Converting to the rest frame of S, 
we obtain the total positron energy spectrum 

x 5 \2x s e x^ F - x 1 / - cose I>F yJ (xl' F Y - Aelsjl - 4^J , (B12) 

with the scaled energies xf = 2E^/rriD and xi t F = ^Ei^/mD, and the mass ratios 
e e = m e /mj and €i t F = mi^/mo. The branch ratio B^p is given by Bjp = Ti^f/T total- No- 
tice that X e denotes the largest positron energy fraction in the fragmentation of the (7, F) 
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particle. But the positron production is not clear in the course of fragmentation, so we 



cannot make sure the value of X e . However, for m e <C rrijp, we have X e = 1 44]. 

For three-body decays, the upper limit Rfp* can be found in Eq. (1B2|) . In the rest frame of 
the (I, F) particle, 9j t p is the angle between the electron spatial momentum direction and the 
boost axis of the (/, F) particle. The 5— function simply indicates that for some fixed boosted 
energy xf, one should integrate over all possible configurations for (Xg' F , cos 9i,f)- In other 
words, Eq. (IB12|) makes the Lorentze boost of the energy spectrum: dN/dx 1 ^ — > dN/dx^, 
and then convoluted by the distribution dN/dxi t p that denotes the differential probability 
of S decay to (J, F) particle in the rest frame of S. Finally, we get the total positron energy 
spectrum. 

For isotropic decay of the (I, F) particle, the distribution function is independent on 
angle, which can be integrated out 

dN T eff Cj F R IiF r R T,r l dN r™*M 



E /=== A dxi ?—„—J., . ^ 

(B13) 



dxf T tota i R _ 4e 2 ~ 7 2 e /iF -R/,F dx I>F J m in(x ItF ) 

1 diV 



(^• F ) 2 - Ael dx 



I,F ' 



The normalization condition of dN/ dx s e in each decay channel is 

/^—^dxf = multij p , (B14) 
dx% 

where multi^p denotes how many electrons are produced by the fragmentation of each 
(I,F) particle. The integrand region in Eq. (1B13I) is given by mm(xj t p) = 1Z^(xpp) and 
max(xj ) i?) = min(X e ,7?. + ) with 

K±{xi,f) = l( x e ± \J l ~ ^ 2 i,F\/(Xe) 2 ~ 4(e e e 7 , F x 7 , F ) 2 ) eJ 2 F xjy (B15) 

If TZ- > X e , the integral is zero. Then the constrained 1Z±, which depends only on the 
ratio r = xf/xi t p, in turn constrains the integral region for Xj t p. If e^p <C 1 is satisfied, 
then R + « x^/(xiej) 3> 1. So the upper limit always takes X e . For example, let us choose 
(J, F) = fi. But in this paper, this is not always true since the sparticles and Higgs bosons 
are heavy. Consequently there are some values of r which gives TZ+ > X e . For the later case, 
one can show that 7Z± is monotonically increase with r, and the lower limit smaller than 
X e gives the upper bound for r < 772, which is the larger root of 1Z- = 0. However, in the 
former case, 1Z- > X e gives bound rji < r < r/ 2 - in both cases, for fixed xf , we must have 
xi > x s e /rj 2 . 
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2. Analytical Spectra from Three-Body Cascade Decays 



As an application, we use the above formulae to calculate the e ± spectrum functions from 
S — > /j, + fermion + scalar followed by fi three-body decay. Notice that C e e < 1 and 
Vi < r < V2, we obtain 2?7 1 e A1 < xf < K^ ax r\i. The spectrum functions dN/dx^ and dN/dx^ 
are respectively given by 



dA 



2(^) 2 (3-2x£), 



dA 



.r, 



(B16) 



where the approximated spectrum is calculated according to Eq. (IB 5 1) . 

If we are only interested in the positrons with energy higher than 1 GeV, we obtain 
7£_~ xf/x^. Integrating over x^ and x^, we have 



dN 
dx^ 



ftma 

5 3 



dx. 



- 3 



1 - x. 



+ 



3 L 



5 - 5xf + 2 (xf) 2 - 2 (xf) 3 + 4 (xf) 3 log 



2 (5-9 (xf) 2 + 4(xf) ; 



(B17) 

The correction term is complicated but important since without it we cannot have the right 
order. This expression is valid for the energy fraction 2e M ?72 < xf < ri2R™ ax . The DM 
three-body cascade decay produces a harder spectrum than the two-body decay since its 
intermediate particle produces a hard spectrum. This is clearly shown in the muon case. 



Appendix C: Fitting Positron Spectra 

In this Appendix, we present the positron spectra got through the above procedure. 
Similarly, we can calculate the anti-proton spectra, although we will not present the details 
here. The basic method is to simulate the positron and anti-proton spectra in the rest 
frame of the sparticles and Higgs bosons via PYTHIA (each state (/, F) must be simulated 
separately and then summed over). We will not pay much attention to the detail fitting but 
present the results directly. Notation is consistent with the one given above, but we change 
the variable to energy, as in the original form from simulation. 

The spectra are divided into three types: the type of sleptons, the type of Higgs bosons 
and the type of Higgsinos. For sleptons, the spectra from the left- and right-handed sparticles 
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have quite different behaviours, since their weak boson decays are different. The total e* 2 
from fragmentations are described by the following functions 



. dN _ 16. 926-15. 60LE 0183 zt.1.438 

fJ'L ■ , r — exp 



dE 
dN 



£i.438 + 14 90 log ^ + Q.05E) 6 - 1 ] - 8.5 x 10~ 5 E° 



55 



16. 913-15. 601.E 0183 7711.437 



£i.437 ^ + 1 49i og ( 1 + 0.055E) 6 ' 5 ] - 9.0 x 10~ 7 E°- 

£1.302 £0.031 



60 



+ 



d£ 0.245 + 262.353£ L9n 104.334 + 5.757 x 10" 5 £ 3 - 242 



(CI) 



The type for all the Higgs bosons have the similar form, probably because they all are 
dominated by the components. For simplicity, we only present A Q as follows 

dN 



A': - = exp 



17.716-15.601£ 0160 7-11.327 



£i.327 j-j + A Q log ^ + Q.021E) 3 - ] - 1.0 x 10~ 6 E - 8 . (C2) 



At last, the fitting functions for the neutralinos C® 4 and charginos Cf 2 are given by 



•<o . 



dN 

~dE 

dN 



Cl, ^1 = exp 15.461-15.601E0 -"3 sl .3ii ^ + ^ Q ^ + Q 02 5£) 2 - 4 ] - 2.7 X lO" 5 ^ 7 , 



C A : — =exp 



: — =exp 



17.880-15.601BO-»4 £1.408 + 14 . 0k)g(1 + 0.024E) 4 ' 2 ] - 1.7 X lO"^ ' 5 , 
17.198-15.601^^ £!.494 ^ _ Wog{l + 0.045£) 3 '°] - 2.7 X lO" 6 ^' 5 , 



Cf: — =exp 



dN 



I7.8i4-i5.60i£°- m E L416 [1 + 13.01og(l + 0.026E) 3 - 5 ] + 1.4 x 10- 6 £ L 



(C3) 



Basically, these fitting functions consist of three parts. The first part describes the low 
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FIG. 6: Typical fitting functions: C® (left), H (middle) and ju^ (right). In the simulations, 3 x 10 5 
events have been generated by PYTHIA. 



energy region very well, the second logarithmic part is added to slow down the exceptional 
decreasing, and the last part is used to modulate the region near the DM mass cutoff region. 
This kind of fitting functions is not universal, and only applies to our parameter sets. We 
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present the typical fittings for three types in Fig. |6j By the way, the small discrepancies 
between the fitting functions and simulations do exist at high energy. So the fitting functions 
are simply used as the referred fitting functions. 
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